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Abstract 

We analyze the phase diagram of QCD with four staggered flavors in the (/j, T ) plane using 
a method recently proposed by us. We explore the region T > 0.7 Tq and p < 1.4Tc, where 
Tq is the transition temperature at zero baryon density, and find a first order transition 
line. Our results are quantitatively compatible with those obtained with the imaginary 
chemical potential approach and the double reweighting method, in the region where these 
approaches are reliable, T > 0.9Tc and p < Tq. But, in addition, our method allows us 
to extend the transition line to lower temperatures and higher chemical potentials. 
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1 Introduction 


The study of QCD at finite baryon density is one elusive problem of utmost importance 
for the understanding of strong interactions. Many features of the behavior of hadronic 
and quark matter at high baryon density have been suggested mainly from the analy¬ 
sis of effective theories. On this basis, a rich structure of phases and phase transitions 
is expected, and the dependence of the phase structure on the quark masses has been 
conjectured |lj. Unfortunately, the sign problem prevents the direct use of Monte Carlo 
simulations at finite baryon density. However, in the last years methods to get the transi¬ 
tion line at high temperature and low baryon density, overcoming the sign problem, have 
been proposed. One method uses double reweighting (in the two parameters, f3 and //a) 
of the configurations generated at the transition temperature at zero chemical potential, 
following the transition line in the (3) plane j2j. The rational of this method is that 
it is expected that the overlap between the reweighted ensemble and the ensemble at a 
given point of the transition line will be improved respect to the Glasgow method J3j, 
since the original ensemble is itself a mixture of configurations corresponding to the con¬ 
fined and deconhned phase. This approach has been used in a first attempt to locate 
the expected critical endpoint of QCD with 2+1 flavors by using only first principles jl|. 
Another method exploits the well known fact that there is no sign problem if the chemical 
potential is purely imaginary 0 . Then, it was realized in UJ that it is possible to deter¬ 
mine, by means of numerical simulations, a pseudo-transition line at imaginary chemical 
potential and extend it analytically to real chemical potential. This method has been used 
to study the phase diagram at small chemical potential around the zero density transition 
point in QCD with two and three degenerate quark flavors in jjjj and with four flavors 
in {7;. Another proposal that is being employed to extract information about the phase 
transition in the region of small chemical potential is to compute the expectation values 
of the derivatives respect to fxa at fia = 0, in order to reconstruct several terms of the 
Taylor series 0Ej. 

Recently we devised another method which can be regarded as a generalization of the 
imaginary chemical potential approach. The former has several advantages over the later, 
for it seems that it can be used to determine the transition line at lower temperatures and 
higher densities m- Especially interesting is the fact that it may be used to locate the 
critical endpoint expected in two flavor QCD En¬ 
in this paper we report the results of the phase transition line of lattice QCD with 
four degenerate flavors of staggered quarks obtained with this new method. At zero 
density a very clear first order transition separates the low and high temperature phases, 
at a transition temperature, Tq, that, for small quark masses, ranges from 100 MeV to 
170 MeV [T2J. It is expected that the transition continues along a line in the (/r, T) plane, 
with the transition temperature lowering as /r increases. This is what has been obtained 
using the double reweighting j2] and the imaginary chemical potential [7] approaches. We 
also find a first order transition line starting at // = 0 and T = Tq and continuing at 
lower temperatures as n increases. Our results are quantitatively compatible with those 
obtained in [7j with the imaginary chemical potential approach and with the results of 
the double reweighting method [2J, in the region where these approaches are reliable (see 
sections three and four for a discussion on this point). But, in addition, our method allows 
us to extend the transition line to lower temperatures and higher chemical potentials. 
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The remaining of the paper is organized as follows: in next section we explain our 
approach. In section three we describe the numerical results and in section four we present 
our conclusions. 

2 Review of the method 

The numerical method used in this work is based on the definition of a generalized QCD 
action which depends on two free parameters (x,y). This generalized action suffers also 
from the sign problem for real values of y but not for imaginary values of it. Simulations 
will be then performed at imaginary values of y and at the end analytical extensions will 
be needed. The main advantage of this approach when compared with the imaginary 
chemical potential method is that we can explore the phase transition line at imaginary 
values of y at any given physical temperature i.e., we are not forced, as in the case of 
imaginary chemical potential, to perform simulations at so high temperatures that the 
system is in the quark-gluon plasma phase for any real value of ya. In this section we 
shall describe this method. The interested reader can find more details of it in m- 
The lattice action for QCD with staggered fermions and chemical potential y is 

1 3 

S = ^pg + ^EE IpnJliijl) (^Un,i^Pn+i 

n i= 1 

T ~ 'y ] '^Pnyoixi) U n y'll ) n +o e M U^ l _Q Q'll> n —Q S j + ma ^ ^ i (1) 

n n 

where Spg is the standard Wilson action for the gluonic fields, which contains j3, the 
inverse gauge coupling, as a parameter, rn{n) and r/o(n) are the Kogut-Susskind phases, 
m the fermion mass, and a the lattice spacing. 

Let us now define the following generalized action 

1 3 - 

S = Spg + mo E V’nV’n + ^ EE tpnViiP’') i^n,i^n+i U n _^ T S T (x, y) , (2) 

n n i= 1 

with 

S T (x, y) = E] ''PnVoin) (Un, oVVi +0 - C4_ 0 ,oV>n—o) 

n 

+ y\ ^2 PnVoin) (u nfi i/j n+ o + c4_o,o^n-o) , (3) 

n 

where x and y are two independent parameters. The QCD action is recovered by setting 
x = cosh (ya) and y = sinh (ya). 

Monte Carlo simulations of the model © for real values of x, y are not feasible since we 
meet the sign problem. However if y is a pure imaginary number, y = iy, where y is real, 
the sign problem disappears since the fermionic matrix is the sum of a constant diagonal 
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Figure 1: Minimal phase diagram conjectured for the generalized QCD. The solid line is line of phase 
transitions. The discontinuous line is the physical line, x 2 — y =1. 


matrix plus an antihermitian matrix which anticommutes with the staggered version of the 
75 Dirac matrix, and numerical simulations become feasible. Imaginary chemical potential 
is a particular case of this, obtained by setting x = cos (pa) and y = sin (pa). 

The expected phase diagram for this model in the (x, y ) plane is shown in Fig. jT] 
HOI. The solid line is a presumed line of phase transitions and the discontinuous line is 
the physical line, x 2 — y 2 = 1, along which one recovers standard QCD at finite baryon 
density. The intersection of the solid line with the discontinuous one will therefore give 
us the transition chemical potential of QCD at a given temperature. A change in the 
physical temperature can be simulated by changing f3 keeping fixed Lt or vice-versa. In 
both cases the solid line in Fig. [Tj will move and the intersection point which gives the 
transition chemical potential will change with the physical temperature. For small (3 the 
transition line crosses the y = 0 axis at x > 1 and, therefore, intersects the physical line 
also at x > 1, producing a physical phase transition at pa > 0. By increasing (3 and 
keeping fixed the temporal lattice extent Lt, the transition point on the y = 0 axis moves 
toward x = 1 and eventually crosses it. Clearly, the value of [3 at which the transition line 
intersects the physical line at x = 1 and y = 0 is the zero density transition point. For 
larger [3 the transition line and the physical line do not intersect and the system is in the 
deconhned phase whatever the chemical potential. 

^From an analysis of the symmetries of the action m it is not difficult to realize that 
the partition function depends on x and y only through the combinations 

u = x 2 - y 2 , 

v = (x T y) 3Lt + (x - y) 3Lt . (4) 

For imaginary values of y (y = iy ) we have u = p 2 and v = 2p 3Lt cos(3L t r]), where 
p 2 = x 2 + y 2 and tan rj = y/x, and, therefore, the free energy will be a periodic function 
of r] with period 2ir/3Lt. In particular if the phase transition line of Fig. 1 continues to 
imaginary values of y, the expected phase diagram in the (x, y) plane is displayed in Fig. 
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Figure 2: Conjectured phase diagram in the ( x, y) plane. We have incorporated the periodicity. The 
dashed line contains the points corresponding to imaginary chemical potential. 


2, where we have incorporated the property of periodicity. In Fig. 2 we have also included 
the line p = 1 (dotted) which is the locus of the points accessible to numerical simulations 
of QCD at imaginary chemical potential. One can see now how this approach has more 
potentialities than the imaginary chemical potential approach. Indeed by increasing the 
inverse gauge coupling (3, the phase transition line of Fig. 2 moves toward the origin of 
coordinates. In some interval (/ 3 m ,PM ) the transition line intersects the p = 1 line and 
then a phase transition will appear at imaginary chemical potential. In such a situation, 
the physical temperature is so high that the system is in an unconfined phase for any real 
value of the chemical potential. The advantage of our approach is that in our simulations 
p is not enforced to be one. 

The variables u and v have the interesting property that are real for both y real and 
pure imaginary, so that they map the two planes (x, y) and (x, y) onto a single plane ( u , v). 
The line v = 2u 3Lt ^ 2 separates the regions corresponding to each plane: the region above 
it in Fig. 3 corresponds to real y and is not accessible to numerical simulations, while the 
region below it corresponds to imaginary y and can be explored by means of simulations. 
The physical line is u = 1 with v > 2. The imaginary chemical potential points are 
mapped onto the line u = 1, —2 < v < 2. The analytical extension from imaginary y to 
the physical real y becomes in the ( u , v ) plane the extrapolation from the region accessible 
to numerical simulations to u = 1 , v > 2 . 

In order to get the transition chemical potential we must determine the coordinates 
of the intersection point between the solid line and the physical line in Fig. |T] ^Frorn the 
symmetries of the partition function we know that the phase transition line is an even 
function of x and y. Therefore, we can write the following equation for the transition line 
at fixed /3 in the (x, y) plane 1 

x 2 = 1 + o 0 (/3) + a 2 (/ 3 ) y 2 + a 4 (/?) y 4 + 0(y 6 ), (5) 

By fixing the lattice temporal extent L t and the gauge coupling f3 one fixes the physical 

1 Here and in the following we do not write explicitly the dependence of the coefficients a,i,j3o,bi , etc. on 
the temporal lattice extent, Lt, and the quark masses, ma. 
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temperature T. The intersection point of the transition line with the physical line 


vl = «o (P) + a 2 (/?) y 2 c + a 4 (/3) y 4 + 0{y%), ( 6 ) 

will give us the transition chemical potential, y c = sinh(// c a), at this temperature. 

The strategy for the determination of the transition chemical potential is then the 
following. From numerical simulations at imaginary values of y, y = iy, near the phase 
transition point ((1 + ao) 1 / 2 , 0 ) one can locate several phase transition points in the (x, y) 
plane (see Fig. 2). By fitting these points with equation © with the + sign of the coeffi¬ 
cient proportional to y 2 replaced by —, we can numerically measure the first coefficients. 
Ignoring the quartic term, the transition chemical potential, y, c a, will then be given by 

y c a = TsinlT 1 ^ ^ ^ ^ . (7) 

An alternative procedure, which is indeed the one employed in this work, is to project 
of the phase diagram onto the (y,/3) plane. In practice, we fix x = xq > 1, L t and the 
lattice quark masses 2 and perform simulations for different values of y and (3. In this way 
we can easily find an accurate estimate of the transition point at a given y by interpolation 
in /3 via reweighting. The qualitative phase diagram in the (y, (3) plane is displayed in 
Fig. 4. The solid lines are the physical lines, 


y = 2/ph = 


± 



( 8 ) 


and the discontinuous line is a line of phase transitions. 

Using again the symmetries of the partition function we can write for this line the 
following equation 


(3 = f3 0 (x 0 ) + b 2 ( x 0 ) y 2 + 0 (y 4 ) . (9) 

As in the case discussed previously, one can measure (3q and b 2 from simulations at 
y = 0 and at imaginary y = iy (keeping xq fixed). Then, we have to find the intersection 
point between the physical and phase transition lines of Fig. 4. 

Equivalently, one may use the plane (u,/3), where u = x^ + y 1 > 1. ^From simulations 
at fixed x = xq and y > 0 we can get a phase transition line (3{u ) for u > Xq > 1. Then 
one can extrapolate this line to u = 1, thus obtaining the physical transition coupling at 
chemical potential /ia = cosh _1 (xo). 

3 Numerical results 

We performed simulations of QCD with four degenerate flavors of staggered quarks using 
the Hybrid Monte Carlo algorithm on lattices of sizes 6 3 x 4 and 8 3 x 4. Our method is 
computationally very expensive and it is rather difficult to go to larger lattices. The quark 
mass in lattice units was fixed to ma = 0.05. On each lattice, we repeated the simulations 
for several values of x , y , and (3. We measured the plaquette, the chiral condensate, 
and the Polyakov loop after each molecular dynamics trajectory of unit time. For each 
simulation we accumulate between thirty and forty thousand measurements. Most of the 

2 Notice that in the imaginary chemical potential approach xo is enforced to be less than one. 
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Figure 3: Phase diagram in the ( u , v ) plane. The solid line, v = 2 u 3Lt ^ 2 corresponds to y = 0, and separates 
the region where numerical simulations are feasible (below the line) from the region where the sign problem 
prevents numerical simulations (above the line). The discontinuous line is a hypothetical phase transition 
line. The line u = 1 is also displayed. For v > 2 it is the physical line while for — 2 < v < 2 it corresponds 
to imaginary chemical potential. 



Figure 4: Conjectured phase diagram in the (y, /?) plane, with x = xo > 1 fixed. The discontinuous line is 
a line of phase transitions, whereas the solid lines are the physical lines, y = y p h = ± 1 / Xq — 1. 


simulations have been performed on the Linux clusters of LNGS-INFN, and some of them 
in the Linux cluster of Departamento de Ffsica Teorica of Universidad de Zaragoza. 

At zero chemical potential (x = 1 and y = 0) there is a very clear signal of a first 
order phase transition controlled by (3 , at which the plaquette, the chiral condensate, and 
the Polyakov loop vary abruptly and show a clear two state structure. This first order 
transition persists for x > 1 and y > 0, with the transition coupling depending on x and 
y. To determine the transition lines numerically, we proceeded as follows. For each pair 
(x, y) we estimate the transition coupling by reweighting a la Ferrenberg-Swendsen (in the 
parameter (3 ) the configurations corresponding to the value of (3 at which a clear two state 
signal appeared. We used the maximum of the plaquette susceptibility as the criterion to 
define the transition coupling on the finite lattice, since it gives the best signal. 

In this way, for each value of x we have a transition line in the plane ( y, (3 ). We fit this 
line with a second order polynomial: 

P(y) = Po(x) + b 2 {x)y 2 . (10) 
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L 

x 

00 

Pph 

b 2 

X 2 /ndof 

ndof 

6 

1.02 

5.0359(2) 

5.008(1) 

0 . 68 ( 2 ) 

0.4 

1 

6 

1.045 

5.028(1) 

4.962(3) 

0.72(3) 

1.63 

1 

8 

1.0201 

5.0367(5) 

5.0089(10) 

0.684(15) 

0.82 

4 

8 

1.0314 

5.0332(3) 

4.9869(12) 

0.726(13) 

1.11 

5 

8 

1.0453 

5.0292(5) 

4.9636(21) 

0.708(21) 

0.97 

5 

8 

1.0811 

5.0179(4) 

4.8943(20) 

0.732(10) 

0.36 

7 

8 

1.1276 

5.0047(4) 

4.8017(28) 

0.747(9) 

0.97 
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Table 1: Parameters of the x 2 fits of the transition lines in the (y,0) [or (u,@)\ plane for 
different values of x, extracted from simulations on the L = 6 and L = 8 lattices. 

The analytic continuation of this functions to real y = —iy is trivial: 

P(y) = 0o(x) - b 2 (x)y 2 . ( 11 ) 

The physical value of y is y p h = \J x 2 — 1. Hence, the physical transition coupling is 
P(y P h) = Po(x) - b 2 (x){x 2 - 1 ). 

We can also use the variable u = x 2 + y 2 , and then the transition line in the plane 
(u, (3) is fit by the linear function 

(3{u) = f3 ph {x) + b 2 (x){u- 1). (12) 

Notice that the fits in the (y,0) and ( u,{3 ) are identical since both the data and the 
fit functions are related by the same transformation. In the (u, (3) plane the analytic 
continuation to y = —iy p h corresponds to the extrapolation to u = 1. Hence, the physical 
transition coupling for each value of x is directly given by the parameter (3 P h. 

Table I collects the parameters of the best x 2 fits of the transition lines given by 
equations m or m extracted from simulations on the 6 3 x 4 and 8 3 x 4. The transition 
lines at fixed x are displayed in Figures El and El The left panels represent the (y, (3) plane, 
for the open symbols and the dashed line, and the (y, (3) plane, for the solid and dotted 
lines, and the filled circle. The right panels represent the (u, (3) plane. In all plots the 
open symbols are the numerical estimates of the phase transition points obtained form 
the simulations on the 8 3 x 4 lattice. On the left panels, the dashed line is the best y 2 fit 
to a function of the form m and the solid line is its analytical continuation to y = iy. 
The dotted line is the physical line, given by y = y p h = Vx 2 — 1, and the filled circle is 
the physical transition point. On the right panels, the solid line is the fit of the transition 
points to a function of the form m and the dashed line is the boundary of the region 
where numerical simulations without sign problem can be performed (i.e., the line u = x 2 ). 
The filled circle is the physical transition point obtained extrapolating the transition line 
to the physical region, u = 1 . 

In the transition lines displayed in Fig. El the points at the largest values of y (or 
u) deviate from the smooth growing of the points at smaller y (or u), and they are not 
included on the fits. They are a reflection of the 2n/3Lt periodicity that we would see 
if we look at the phase diagram onto the constant f3 plane instead of onto the plane of 
constant x (see section 2). The data of Fig. 6 (especially those for x = 1.0811, for which 
many points of the phase transition line were determined) strongly suggest that there is 







a cusp separating the /3 increasing curve from the /? decreasing one. This means that 
the 27 t/3 Lt periodicity is realized by the periodic replication of a nonperiodic analytic 
function, producing a cusp at the points where a replica ends and a new replica starts. 
This cusp is a non-analyticity of the phase transition line which has the same origin as 
the Roberge-Weiss transition m found in QCD at y = ±i7r/3Li in the high temperature 
region M- Indeed, as it was conjectured in El, this singular behavior realized as a 
cusp is to be expected in a wide variety of models characterized by a quantized charge 
coupled to a phase. In our case the quantized charge is the baryonic charge and the phase 
(f) = arctan(y/x). Under this conditions, one must exclude the points to the right of the 
cusp from the fits since they belong to a different analytic function. Notice also that the 
presence of the nonanalytic cusp does not necessarily imply that the convergence radius 
of the Taylor series around y = 0 (u = 1) is limited by this singularity. Reference IJj 
contains a simple illustrative example on that in what we called gaussian model. 

The transition lines displayed in Figures El and El are very smooth [they are essentially 
straight lines in the (u, (5) plane until the cusp] and suggests a linear fit in the (u, /3) plane. 
In fact the data reported in Table I imply actually a very high confidence level for the fits 
and to add higher order corrections, as for instance a quartic term in equation © , seems 
meaningless. In other words, if higher order corrections were relevant at real y, it would 
be very hard to measure them from the data produced at imaginary y. This makes really 
difficult any serious analysis of systematic errors possibly induced by the fit ansatz. 

Figure El displays the phase diagram in the plane (ya,/3), together with the results 
obtained by Fodor and Katz with double reweighting j2j, and by D’Elia and Lombardo 
with imaginary chemical potential simulations jjj. There is good agreement with the 
results of Fodor and Katz until ya ~ 0.3, especially if we take into account that different 
methods to locate the pseudotransition coupling at finite volume are used. The agreement 
with D’Elia and Lombardo is very satisfactory in the whole range of ya explored. The 
main difference, which is already seen at ya = 0, should be attributed to a volume effect, 
since these authors used a 16 3 x 4 lattice. D’Elia and Lombardo give credit to their 
results for ya < 0.3, which is the interval where they found agreement with Fodor and 
Katz. We have seen that indeed the results (ours and theirs) are reliable at least until 
ya = 0.5, which is the maximum chemical potential at which we performed simulations. 
The disagreement with Fodor and Katz for ya > 0.3 is likely due to the poor overlap of 
the ensemble of reweighted configurations with the ensemble of typical configurations at 
ya > 0.3. 

Figure HO shows the phase diagram in the plane (y,T) in physical units, with the scale 
set by the transition temperature at y = 0, Tq. The relative lattice spacings have been 
determined by means of the two loop beta function. For comparison, we also plot the 
result of D’Elia and Lombardo. Our results can be fit with a power function of the form 

= [1 -c{y/T c f] p . (13) 

The best \' 2 fit gives c ~ 0.446(9) and an exponent, p ~ 0.173(7). One may be curious 
about the extrapolation of this line to zero temperature. If this were done, we would find 
a zero temperature transition at chemical potential yc ~ 1.5Tc, which, in terms of the 
nucleon mass El, tun, is yc ~ tttn/5. 
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4 Conclusions 


Using a recently proposed method to determine the phase transition lines at finite baryon 
density, we have obtained the phase transition line of QCD with four degenerate quark 
flavors from the zero density high temperature transition, at Tq, down to T ~ 0.7Tc 
and g, ~ 1.4Tc. Our results are in reasonably good agreement with those of D’Elia and 
Lombardo, obtained from simulations at imaginary chemical potential. These authors give 
credit to their results for T > 0.9Tc, since this is the region where their results showed 
reasonably small statistical and systematic errors and in addition agree with those obtained 
by Fodor and Katz with the double reweighting method. Our results agree reasonably well 
with the central value reported by D’Elia and Lombardo in the whole region T > 0.7Tc 
and n < IAT C . 

We believe we can explore much lower temperatures with our method. Since it is 
based on analytical continuation/extrapolation there are uncontrolled systematic errors 
that grow in decreasing the temperature. It is therefore very difficult to estimate the 
minimum temperature at which our method will give a reliable prediction of the location 
of the phase transition point. In m we verified that in the three dimensional Gross-Neveu 
model at large N this minimum temperature is amazingly low. 

The study of four flavor QCD thermodynamics at lower temperatures is very interest¬ 
ing, and we left it for future work. Since our method is computationally very expensive, 
we decided to concentrate the present effort in the more interesting cases of two and two 
plus one flavor QCD. 
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Figure 5: Phase diagrams at fixed x in the (y, /3) 
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Figure 6: Phase diagrams at fixed x in the (y, j3) and (y,/3) planes (left panels) and in the (u,/3) plane 
(right panels). 
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Figure 7: Phase diagram in the (fia, (3) plane. The solid line is the analytical continuation of the imaginary 
chemical potential pseudotransition line, and the dashed lines mark the error band. The authors of Ref. JJ 
give credit to this analytical continuation up to [ia « 0.3. 
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Figure 8: Phase diagram in physical units, with the zero density transition temperature, Tq, setting the 
scale. 
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